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METHOD AND APPARATUS FOR REDUCING IMPULSE NOISE IN A SIGNAL 

PROCESSING SYSTEM 

BACKGROUND OF THE INVENTION 

1. Field of the Invention 

The present invention relates generally to signal processing, and more particularly, to 
a method and apparatus for reducing impulse noise in signals transmitted using 
communication services or recorded using imaging devices. 

2. Description of Related Art 

Currently, there is a significant desire to exploit the unused available bandwidth of the 
twisted pair lines of the existing plain old telephone system (POTS) for providing various 
digital services. Although it is believed that the future media for networked data transmission 
will be fiber optic based and although the main backbone of the network that interconnects the 
switching centers is now mainly optical fiber, the 'last mile' which is the access portion of the 
network that connects switches to customers is still dominated by twisted copper wires. For 
example, there exits over 560 million 'last mile' twisted copper pair connections globally. The 
estimated cost of replacing these connections with fiber optics is prohibitive and therefore the 
existing unused bandwidth of the POTS provides an important alternative. 

Advanced digital transmission techniques such as digital subscriber line services 
utilize the existing unused bandwidth of the POTS for providing increased data transmission 
rates for available digital data transmission services. By convention, "digital subscriber line" 
services are referred to as "DSL" services. The term "DSL" refers a connection created by a 
modem pair enabling high-speed digital communications. More generally, DSL is referred to 
as xDSL, where the 'x' indicates a number of different variants of the service (e.g., H (High), S 
(Single-Line), and A (Asymmetric)). 

One factor that impairs the performance of xDSL services or other similar services 
that operate at high frequencies, such as integrated digital services network (ISDN), is 
"impulse noise." Impulse noise is noise that occurs with high amplitudes on telephone lines or 
other transmission mediums. That is, samples of impulse noise have very large amplitudes 
that occur much more frequently than they would with Gaussian noise. Some known causes 
of impulse noise include electrical equipment operating near the telephone line or relay re- 
openings and the ringing of a telephone on the line. 

In operation,' xDSL services rely on modems to carry digital signals over the pass- 
band channels of the POTS. The modems translate digital data to analog signals at the 
sender end of the telephone line and translate the analog signals to digital data at the receiver 
end of the telephone line. The analog signal output at' the receiver end of a telephone line is a 
corrupted version of the analog signal input at the sender-end of the telephone line. 
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More specifically, the analog signal output from a telephone line is generally referred 
to as an "observed" signal. The observed signal includes a noise component and data 
component. An observed signal without the noise component is defined herein as a clean 
signal. In order to recover the data component from the observed signal, impulse noise 
introduced during the transmission of the data component must be identified. 

One technique for recovering the data component is to estimate (i.e., predict) what the 
clean signal is without the noise component. Data components of output signals that are 
estimated are referred to herein as "cleaned" signals. One such estimation technique isolates 
the noise component from the data component in an observed signal by modeling the noise 
component using a probability density function (i.e., pdf) that describes the observed 
statistical properties of the noise component. 

Once the noise component is accurately modeled using a pdf, the pdf can be used to 
define an error criterion (also referred to herein as a cost function). The error criterion is 
minimized to solve for model parameters, which are used to estimate the data component of a 
sampled signal. 

A common pdf used to model noise is a Gaussian (or normal) distribution. One factor 
for using a Gaussian distribution to estimate noise is that the Gaussian assumption leads to 
simple estimation techniques. The reason the Gaussian distribution does not accurately 
estimate impulse noise is because impulse noise exhibits large amplitudes known as outliers 
that occur too frequently to fit to a Gaussian model. This characteristic suggests that the 
underlying probability distribution that models the noise has heavier tails as compared to a 
Gaussian distribution. 

It has been suggested that an alpha-stable distribution is one alternative to a 
Gaussian distribution for modeling impulse noise. Because there exists no compact form to 
express its probability distribution function, an alpha-stable distribution is typically defined by 
its characteristic function cp(z), which is the Fourier transform of its probability density function. 

<p(-) = exp{j&~ r \z\ a [\ + jj3s\gn(i)w(z,a)]} (1) 
where, 

a is the characteristic exponent such that 0 < a < 2, 
p is the symmetry parameter such that -1 < (3 < 1 , 
y the dispersion such that y > 0, 
5 is the location parameter such that -<x> < 6 < oo, and 





tan , if a * I 



vt<r,a) = < 



hog\z\, ifa = \ 



More specifically, the parameters control the properties of the pdf of an alpha-stable 
distribution as follows. The characteristic exponent a is a measure of the thickness of the tails 
of the alpha-stable distribution. The special case of a=2 corresponds to the Gaussian 
distribution, and the special case of a=1 with (3=0 corresponds to the Cauchy distribution. The 
symmetry parameter p sets the skewness of the alpha-stable distribution. When p=0 the 
distribution is symmetric around the location parameter 5, in which case the alpha-stable 
distribution is called a symmetric alpha-stable (i.e., SaS) distribution. The location parameter 
5 determines the shift of the alpha-stable distribution from the origin, and is the mean (if 
1<a<2) or median (if p=0) of the alpha-stable distribution. Finally, the dispersion y measures 
the deviation around the mean in a manner similar to the variance of a Gaussian distribution. 

Alpha-stable distributions have been used to design systems for detecting signals in 
the presence of impulse noise. (See for example, E.E. Kuruoglu, W.J. Fitzgerald and P J.W. 
Rayner, "Near Optimal Detection of Signals in Impulsive Noise Modeled with a Symmetric 
alpha-Stable Distribution", IEEE Communications Letters, Vol. 2, No. 10, pp. 282-284, 
October 1998.) However, most of these systems that use alpha-stable distributions in their 
statistical models, assume a priori knowledge of the parameters of the alpha-stable 
distribution. Systems that assume a priori knowledge of the parameters of an alpha-stable 
distribution pre-assign values for the parameters. Having the ability to estimate, and not pre- 
assign, the value of parameters of the alpha-stable distribution is vital since most existing 
systems are sensitive to the parameters of the alpha-stable distribution that models the 
impulse noise. 

Existing methods for estimating parameters of an alpha-stable distribution generally 
provide limited solutions for the special case of a symmetric alpha-stable distribution (SaS) 
(i.e., where the parameter p=0). Assuming that an alpha-stable distribution is symmetric, 
however, may yield a poor model of impulse noise because impulse noise tends to be more 
accurately modeled by skewed rather than symmetric distributions. Existing methods for 
estimating the parameters of an alpha-stable distribution, which provide general solutions that 
are not limited to the special case of a symmetric distribution, tend to be computationally 
expensive or provide estimates with high variances. 

It would be advantageous therefore to provide an improved system for modeling 
additive impulse noise corrupting data streams. Furthermore, it would be advantageous if 
such a system were able to model impulse noise using an alpha-stable distribution. Also, it 
would be advantageous if the improved system were able to adaptively estimate, and not pre- 
assign, the parameters of an alpha-stable distribution. 
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SUMMARY OF THE INVENTION 

In accordance with the invention, there is provided a signal processing system for 
reducing impulse noise corrupting sampled signals. A memory of the signal processing 
system accumulates sampled signals from a transmission medium. The sampled signals have 
a noise component and a data component. In one embodiment of the invention, signals are 
sampled after being transmitted over a transmission medium such as a digital subscriber line 
(DSL) service. In another embodiment of the invention, signals are sampled from a 
transmission medium such as a sensor array in an imaging system such as a scanner. 

In accordance with one aspect of the invention, a parameter estimation module 
estimates the parameters of an alpha-stable distribution. The alpha-stable distribution is used 
to model impulse noise corrupting data signals input into the transmission medium of the 
signal processing system. A coefficient optimization module uses a modified iteratively 
reweighted least squares (IRLS) technique to optimize the model coefficients of a prediction 
filter, such as a Volterra filter. Using the model coefficients, the prediction filter computes an 
estimate of the data component of the signals sampled from the transmission medium without 
the noise component. 

In accordance with another aspect of the invention, the parameters of an alpha-stable 
distribution are estimated using a sampled signal having only a noise component. In the 
embodiment in which the signal processing systems operates a DSL service, a clean signal is 
transmitted over an analog data channel. To sample a signal without a data component, the 
analog data channel is sampled when no data signals are transmitted over the data channel. 
In contrast, in the embodiment in which the signal processing system operates in an imaging 
system, a sampled signal containing only a noise component is generated by applying centro- 
symmetrizing and centralizing transformations to corresponding pixels from multiple recorded 
images of the same scene. 

In accordance with yet another aspect of the invention, the characteristic exponent of 
an alpha-stable distribution is used to define the order of the moment in the cost function that 
optimizes estimation of cleaned signals by the prediction filter. In effect, the cost function is 
defined to be the p th -power error criterion, and the modified IRLS technique is applied to 
optimize the model coefficients of the prediction filter. 

Advantageously, the present invention provides a method and apparatus therefor, for 
modeling impulse noise in xDSL services using an alpha-stable distribution. In addition, a 
number of different methods for computing parameters of the alpha-stable distribution are 
disclosed. 'Generally, these different methods for estimating parameters of an alpha-stable 
distribution include the steps of performing transformations and computing moments. 



BRIEF DESCRIPTION OF THE DRAWINGS 

These and other aspects of the invention will become apparent from the following 
description read in conjunction with the accompanying drawings wherein the same reference 
numerals have been applied to like parts and in which: 

Figure 1 illustrates an operating environment of a signal transmission system for 
performing the present invention; 

Figure 2 illustrates a detailed block diagram of the noise suppression module shown 
in Figure 1; 

Figure 3 illustrates a general block diagram that represents the different elements 
forming a Volterra filter; 

Figure 4 illustrates an example of a linear filter that forms part of the Volterra filter 
shown in Figure 3; 

Figure 5 illustrates an example of a quadratic filter that forms part of the Volterra filter 
shown in Figure 3; 

Figure 6 illustrates an example of a cubic filter that forms part of the Volterra filter 
shown in Figure 3; 

Figure 7 is a flow diagram that sets forth the steps for adaptively determining the 
coefficients of the Volterra filter using a modified iteratively reweighted least squares (IRLS) 
technique; 

Figure 8 illustrates a flow diagram that sets forth the steps that are performed by the 
parameter estimation module to estimate the parameters of an alpha-stable distribution; 

Figures 9A and 9B are flow diagrams which set forth two different combinations of 
steps for estimating the parameters a, P, y, and 5 of an alpha-stable distribution; 

Figure 10 illustrates an alternate operating environment of a digital image processing 
system for performing the present invention; 

Figure 11 illustrates a process for cleaning impulse noise from digital images in 
accordance with the alternate operating environment of the present invention; and 

Figure 12 is a flow diagram that sets forth the steps performed by the pure noise 
extractor to produce an observed signal block that consists of impulse noise absent of image 
content. 



DETAILED DESCRIPTION 

A. Operating Environment 

Figure 1 illustrates an operating environment of a signal transmission system (i.e., 
signal processing system) for performing the present invention. The operating environment of 
the signal transmission system includes a multi-functional device 102 that communicates with 
other devices over a broadband network 104. The multi-functional device 102 receives and 
transmits digital data over the broadband network through a central office terminating unit 106 
and a remote terminating unit 108. The terminating units 106 and 108 form a modem pair that 
operate together to transmit digital data over an analog data channel (or pass-band channel) 
110. In one embodiment, the analog data channel 110 is a twisted pair line of the plain old 
telephone system (POTS). 

The terminating units 106 and 108 have switches 118. Each of the switches 118 have 
two operating positions A and B. In the operating position A, the terminating units are in 
normal operating mode, during which digital data is transmitted between the multi-functional 
device 102 and the broadband network 104 over the analog data channel 110. In the 
operating position B, the terminating units couple the input to the analog data channel with a 
null modem 116. The purpose of the null modem 116 is to sample the analog data channel 
1 1 0 when it is absent of data signals. As discussed in more detail below, the null modems 1 1 6 
provide the noise suppression modules 114 with a sampled signal consisting of a noise 
component and no data component. 

In accordance with one aspect of the invention, the central office terminating unit 106 
and the remote terminating unit 108 operate together to provide a digital subscriber line 
(xDSL) service. Each of the terminating units 106 and 108 includes a modem 112 for 
transmitting digital signals over the analog data channel 110. The modems 112 receive 
signals filtered by a noise suppression module 114. The noise suppression module reduces 
impulse noise corrupting signals transmitted over the analog data channel 110. To transmit 
and receive digital data, the modems 112 in the terminating units 106 and 108 typically 
include a modulator unit and a demodulator unit. To transmit digital data, the modulator unit of 
a modem receives digital data and encodes the digital data into one or more symbols having a 
plurality of bits. Each encoded symbol is then input to a transmit filter which is used to 
produce a continuous time signal. The continuous time signal is transmitted over the analog 
data channel 110. 

The signal sampled at the output of either end of the analog data channel 110 is 
defined herein as the observed signal x ( . In accordance with another aspect of the invention, 
the observed signal x t is processed by the noise suppression module 114 before being 
demodulated by a demodulator unit and decoded by a decoder unit in the modem 112. The 
demodulator unit of the modem receives a cleaned signal which is the output of the noise 



suppression module 1 14. The symbols output by the demodulator unit of the modem are then 
input to the decoder unit of the modem to produce digital data. When the modem 112 forms 
part of the remote terminating unit 108, the multi-functional device 102 receives the digital 
data output by the decoder unit of the modem. Alternatively, when the modem 112 forms part 
of the central office terminating unit 106, the digital data is output to the broadband network 
104. 

B. Overview Of The Noise Suppression Module 

Figure 2 illustrates a detailed block diagram of the noise suppression module 114 
shown in Figure 1. The noise suppression module 114 corrects distortions caused by impulse 
noise introduced to analog signals propagating along the data channel 110. The 
characteristics of additive impulse noise corrupting input signal x t are typically unknown. 
Consequently, the elements of the noise suppression module 114 are used to estimate (i.e., 
predict) what the output signal x t is without noise (i.e., the cleaned signal Generally, the 
elements of the noise suppression module include a data latch 202, a signal estimation 
module 200 and a parameter estimation module 206, In one embodiment, the signal 
estimation module 200 includes a noise symmetrizer 216, a non-linear prediction filter 210, 
and a coefficient optimization module 208. 

In operation, an observed block of L signals x t is input to the noise suppression 
module 114 and stored in the data latch 202. The signals forming the observed block of 
signals are sampled at some predetermined interval from the analog data channel 110. The 
data latch 202 is a memory which stores L sampled data signals output from the analog data 
channel 110. The signals input to the noise symmetrizer 216 and the non-linear prediction 
filter 210 are delayed by one block of signals (i.e., where each block of signals has a 
length of L samples. The observed block of signals which is stored in the memory of the data 
latch 202 can be represented in a matrix form as follows: 



x[txL] 
x[txL + L-\] 



where t = 0, 1, 2 ; 3, 



The parameter estimation module 206 estimates one or more of the parameters a, (3, 
y; and 5 of an alpha-stable distribution, which are defined above in equation (1). As shown in 
Figure 1 , signals are sampled without a data component by one of the terminating units 106 or 
108 when the switch 118 is set to operating position B. When the switch 118 is set to 
operating position B, the terminating units 106 and 108 are in parameter estimation mode. In 
operating position B, null modems 116 are used to insure that no data signals are output to 
the analog data channel 110 so that an accurate measurement of the noise on the analog 
data channel can be performed. In contrast, when the switch 118 is set to operating position 
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A t the terminating units 106 and 108 are in signal estimation mode where signals received by 
the noise suppression module 114 are used to estimate a clean signal (i.e., ^). 

In one embodiment, a measure of the noise on the analog data channel 110 is taken 
and the parameters a, p, y, and 6 are estimated once and hard-coded or fixed as input to the 
coefficient optimization module 208. In an alternate embodiment, the parameters a, p, y, and 5 
are adaptively estimated and modified during the operation of the noise suppression module 
using a new observed block of signals In this alternate embodiment, the switches 118 
transition from operating position A to operating position B to record samples of noise on the 
analog data channel 110 thereby momentarily interrupting transfer of data traffic transmitted 
over the analog data channel 110. 

Once estimated, the parameters a, p, y, and 5 of the alpha-stable distribution are input 
to the coefficient optimization module 208. In one embodiment, the coefficient optimization 
module 208 optimizes model coefficients a, b t and c of the non-linear prediction filter 210 
using a modified iteratively reweighted least squares (IRLS) technique. The model coefficients 
£, b % and c are then input to the non-linear prediction filter 210 to estimate what the observed 
signal block x t .j is without impulse noise. The cleaned signal block v t% which is an estimate of 
the signal block .vy_/ without impulse noise, is defined in matrix form as: 



v = 
— / 



y[txL] 
>{/ x L + L - 1] 



where t - 0, 1, 2, 3, ... 



In one embodiment, the non-linear prediction filter 210 is a one-dimensional (i.e., 1-D) Volterra 
filter. Those skilled in the art will appreciate that the non-linear Volterra prediction filter 210 
has a non-linear dependence on its input data and a linear dependence on in its coefficients a, 
b y and c. Volterra filters are known in the art as disclosed by M. Schetzen, The Volterra and 
Wiener Theories of Nonlinear Systems, New York: John Wiley & Sons, 1980. It will also be 
appreciated by those skilled in the art that the non-iinear prediction filter operates in an 
extrapolatory mode (i.e., extrapolation). The extrapolatory mode involves the prediction of 
future values using observations from past time steps (i.e., predicting values at time t=T, using 
observations having time steps at time t<T). 

In alternate embodiments of the non-linear prediction filter 210, other non-linear filters 
that are linear in their coefficients such as Radial Basis Function filters (which are known in 
the art as disclosed by B. Mulgrew, in "Applying Radial Basis Functions," IEEE Signal 
Processing Magazine, Vol.13, No.2, pp.50-65 March 1996) and Self-Exciting Threshold 
Autogregressive (SETAR) filters (which are known in the art as disclosed by H.L Koul and A. 
Schick, in "Efficient Estimation In Nonlinear Autoregressive Time-Series Models," Bernoulli, 
1997, Vol.3, No.3, pp.247-277) are used in place of a Volterra filter. 
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In another alternate embodiment of the non-linear prediction filter 210, the non-linear 
prediction filter operates in an interpolatory mode (i.e., interpolation) rather than an 
extrapolatory mode. In the interpolatory mode, observations from both past and future time 
steps [t-k ( t-k+1 l ...,t-1,t+1,t+2...] are used to predict the value of the data at time step t. It will 
be understood by those skilled in the art that this alternate embodiment results in the same 
formulas as presented here in Section C up to a relabeling of the time step indices. For 
example given eight observed signals^c t=0 , * t=1 ; x t=2l ,r ts3 , x t=5 , x t=6 , and .v t=7l the data component 
v tx4 of a signal .y 1s4 is estimated (i.e., predicted) using the eight observed signals. 

In addition to the parameters a, p, y, and 6 of the alpha-stable distribution, the 
observed signal block x ( _/, and the extended matrix X ext are input to the coefficient 
optimization module 208. As described in more detail below, the coefficient optimization 
module 208 uses the parameters of the alpha-stable distribution to specify an l p -norm 
estimation criterion (i.e., cost function). The cost function is minimized by the coefficient 
optimization module 208 to determine the model coefficients a, b, and c of the non-linear 
prediction filter 210. However, because the l p -norm estimator only produces unbiased 
estimates when the noise in the observed signal block x t „i is symmetric, the noise 
symmetrizer 216 may be required to deskew and centralize the noise in an observed signal 
block x t . In an alternative embodiment, a zero th order (i.e., constant) (e.g., qq in equation (2) 
below) term is included in the Volterra filter to compensate for bias in the l p -norm estimation. 

C. Non-Linear Prediction Filter 

The non-linear prediction filter 210 uses model coefficients a, b, and c to estimate the 
cleaned signals v/. The model coefficients a, b x and c are optimized using a parameter of an 
alpha-stable distribution that models impulse noise corrupting the observed signal x t .]. A 
general alpha-stable distribution is different from the Gaussian distribution because the alpha- 
stable distribution lacks finite second order statistics. As a result, the prediction filter 210 
cannot use conventional least squares estimation techniques that are based on minimum 
mean squared error criterion to accurately estimate the cleaned signals v t , since such 
techniques employ second order statistics. 

It is known that minimizing the dispersion of a parameterized random variable 
distributed with an alpha-stable probability density function is equivalent to minimizing the p th 
order moment of the random variable's probability distribution (see for example V.M. 
Zolotarev, "Mellin-Stieltjes Transforms In Probability Theory," Theory of Probability and 
Applications, vol. 2, no. 4, pp. 433-460, 1957). Whereas the minimum mean squared error 
criterion leads to least squares estimation (l 2 -norm), the minimum mean p th -power error 
criterion leads to l p -norm estimation . 

Although the minimum, mean squared .error criterion leads to a linear predictor for 
Gaussian data with Gaussian noise, the error criterion for alpha-stable data or alpha-stable 



noise need not be linear. The filter 210 is, therefore, selected to be a non-linear Volterra filter 
or polynomial filter even if the process generating the clean data may be modeled as linear. 
The non-linear Volterra filter is used to estimate the data component of the observed signal x t _ 
l or of the centro-symmetrized observed signal xj_/. Those skilled in the art will appreciate 
that the observed signal need not be centro-symmetrized before being input to the non-linear 
prediction filter if the noise in the observed signal is symmetric or if zero th order terms are 
included in the Volterra filter. The estimate of the data component of the observed signal 
is defined herein as the estimated cleaned signal yj. 

Using the estimated cleaned signal £/, the noise signal (or component) of the 
observed signal x t can be estimated using an additive model which assumes that the noise 
signal is produced independently of the data signal (or component). The estimate of the noise 
signal is defined herein as the estimated noise signal r t . The relationship between the 
observed signal */_/, the estimated cleaned signal and the estimated noise signal r t can 
therefore be represented using the additive model as: 

-I/-/ = y_t + r j- 

The input-output relationship of a Volterra filter can be defined as: 

N N N 

y(n) = a 0 + Z OjX{n - / ) + Z Z b f , x(n - /) x(n - j) + 

i=\ i=\j=i 

jV jV jV 

Z Z Z c './.* x{n - - - *> + - ( 2 ) 

where: 

x(n) is the observed signal, 

y(n) is the data component or cleaned signal ( 

n is the observed block index that runs from 0 ... L-1, and 

a 0 , a it b ip and c iiM are the Volterra filter coefficients. 

In one embodiment, the Volterra filter 210 is defined for computational efficiency with only the 
first three terms (excluding the term a 0 ) of the general Volterra filter set forth in equation (2). 
Limiting the general Volterra filter to its the first three terms defines a truncated Volterra filter 
having up to third order non-linearity. It will be appreciated, however, by those skilled in the art 
that in alternate embodiments the filter 210 can be defined using truncated Volterra filters that 
have less than or more than three terms. 

Using this input-output relationship, the data signal % for 3 signal block is computed 
given the observed signal block x t .j and the model (or Volterra) coefficients a, 6, and c. The 
model coefficients a % b, and c are received from the coefficient optimization module 208. In 
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computing the data signal the observed signal block x t .j is delayed by one sample. In 
operation, the Volterra filter uses a block of signals sampled at the times [Lxt-L-1+k, Lxt-L+k, 

.... Lxt-2+k] to estimate what a cleaned block of signals is at the times [Lxt-L+k, Lxt+k 

Lxt-1+k] for some k in the range 0 to N-1. The block length L is chosen to be substantially 
longer than the number of coefficients in the Volterra filter. 



form as: 



The input-output relationship of the Volterra filter can also be represented in matrix 
y = X ext C where, 

the extended Volterra data matrix^, = [x 0) X {2) X (3) J, such that: 



x[txL] 



0 

x[( x L] 



x[txL + L~\] x[txL + L-2] x[t*L + L-N] 



x z [txL] 
x 2 [txL + \] 



0 



x[tx L + \]x[txL] 
x 2 [t xL + N] x[i x L + N]x[t xL + N-l] 
x 2 [t xL + L - 1] x [t x L + L - \\x\t xL + L-2] 



,V (3 > = 



x'[txL] 
x 3 [txL + \] 



x~[txL + \]x[txL] 



x 3 [txL + N] x 2 [txL + N]x[t xL + N-\] 



0 
0 

x 2 [txL] 
x 2 [txL + L-N] 

0 
0 

x 3 [txL] 



_x 3 [t x L + L - 1] x 2 [txL + L- \]x[t xL + L-2] * - * x 3 [txL + L-N] 
where: 

each row in x( 2 ) (and x(V) corresponds to the quadratic (cubic) terms in the Volterra 
expansion for a fixed time instant t given in equation (2) above, 

t is a fixed time constant having a range from 0,1, ... , N, where N represents the 
memory capacity of the filter 210, 

L is the data block length size; and 



the Volterra vector of coefficients C = 



, such that: 
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a \ 








^1 


a = 




, b = 


2 


, and c = 


6 



Figures 3-6 illustrate different elements of a non-linear Volterra filter 210. Figure 3 
illustrates a general block diagram that represents the different elements forming a Volterra 
filter/ As illustrated, the delayed observed signal block x t .j is input to a linear filter 302, a 
quadratic filter 304, and a cubic filter 306. The output of the filters 302, 304, and 306 is the 
product of the Volterra data matrices and Volterra coefficients defined above as xfVa, x(2)b, 
and x(*)c. The elements in each of the resulting vectors are summed together by summation 
unit 308 to provide the estimated cleaned data signal yj. 

Figure 4 illustrates one embodiment of a linear filter 302. In operation, the linear filter 
302 has shifted through the first register of shift register 402 each element in the sequence (or 
vector) of observed data signals starting with the first element of the block x t .j{n=0). After 
each shift of the shift register 402, the coefficients are multiplied by the Volterra a coefficients 
using multipliers 404. These results are summed at adders 406 and output to define the 
entries in the resulting vector x(V a. Figure 3 illustrates the computation of the first element in 
the vector JtfOa. This operation can be summarized as the convolution of the input signal 
block ,y/_/ and the impulse response of the filter 210 (i.e., the filter coefficients a, b t and c). In 
an alternate embodiment, the operation in Figure 3 is performed using overlapping blocks. 

Figure 5 illustrates an example of the quadratic filter 304 illustrated in Figure 3. The 
quadratic filter 304 can be defined using a quadratic sequence generator 502 and a linear 
filter 504. The output of the quadratic sequence generator is a quadratic sequence which is 
input to the linear filter 504 to produce the resulting vector x( 2 )b which can be viewed as the 
convolution of the quadratic sequence with the filter coefficients b. Similarly, Figure 6 
illustrates an example of a cubic filter 306 illustrated in Figure 3. The cubic filter 306 includes 
a cubic sequence generator 602 for generating a cubic sequence. The resulting cubic 
sequence is subsequently input to the linear filter 604 to produce the resulting vector xO)c. 
Each sample in the quadratic (cubic) sequence corresponds to (or generated such that they 
correspond to) a quadratic (cubic) term in the Volterra filter expression given above in 
equation (2). One reason for representing the Volterra filter in matrix form is to simplify the 
computations of the coefficient optimization module 208 by utilizing the extended matrix X ext - 

D. Coefficient Optimization Module 

As illustrated in Figure 2, the coefficient optimization module 208 receives as input 
parameters from the parameter estimation module 206 and the extended matrix X ext from the 
non-linear prediction filter 210. The coefficient optimization module 208 adaptively determines 
the coefficients of the Volterra filter using a modified iteratively reweighted least squares 
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(IRLS) technique. The steps for performing this technique are summarized in the flow diagram 
set forth in Figure 7. 

Initially at step 700, the index k is initialized to zero. Also at step 700, the weight 
matrix W is initialized to an identity matrix /, and the value of \\rf-I)\\ {p) is initialized to zero. At 
step 702, the value of p is set equal to the value of the characteristic exponent a received 
from the parameter estimation module 206. In accordance with this aspect of the invention, 
the value of the characteristic exponent a is used to define the order of the moment used to 
compute the model coefficients of the non-linear prediction filter 210. 

At step 704, an initial value for the vector of Volterra coefficients C(0) is computed for 
k=0. Subsequently, at step 706, an error signal r/ (or residual error term) is computed for each 
/ in (0...L-1) using the observed signal block xt-I> the extended Volterra data matrix X ex i and 
the vector of Volterra coefficients C(k). At step 708, elements Wa of the diagonal weight matrix 
W are computed for each / in (0...L-1). The resulting vector of error signals r t and diagonal 
weight matrix Ware defined as: 



and W = 



o,o 
0 

0 



0 



0 



0 IV i 



£-1,1-1 



At step 710, a vector of Volterra coefficients C(k+1) is computed for the subsequent 
index value (e.g., k+1) using the computed diagonal weight matrix, the extended Volterra data- 
matrix X exh and the observed signal block x t . At step 712, a determination is made as to 
whether the error criterion for estimating the Volterra coefficients has sufficiently converged. 
Sufficient convergence is achieved when the relative change in the norm of the estimation 
error \\tj\ {p) between iterations is smaller than the convergence limit s. In one embodiment, the 
convergence limit e equals 10" 4 . The error criterion \\r\\ {pY which is the p th -power error criterion, 
is computed as follows: 

When convergence is successfully achieved, then step 716 is performed and the 
vector of Volterra coefficients C(k+I) last computed at step 710 is passed to the non-linear 
prediction filter 210. If the solution did not successfully converge then step 714 is performed. 
At step 714, the index k is incremented and step 706 is- repeated. It will be appreciated by 
those skilled in the art that an upper limit of the index k can be defined in order to assure that 
a Volterra coefficient vector is found at step 716 within a constrained amount of time. 

E. Noise Symmetrizer 
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In general, the coefficient optimization module 208 can only produce unbiased 
estimates of the coefficients of the Volterra model with no zero 1 * 1 order term if the impulse 
noise has a symmetric probability density function. In accordance with another aspect of the 
invention, the noise symmetrizer 216, which includes a random noise sequence generator 212 
and a differencer 214, is adapted to convert observed signal blocks with impulse noise having 
non-symmetric probability density functions into a form that can be used to compute an 
unbiased estimate of the coefficients of the non-linear prediction filter 210. This aspect of the 
invention relies on the assumption that there exists a means for obtaining replicas of observed 
signal blocks with the same data component but different noise component that are derived 
from the same statistical distribution. 

More specifically, the random noise sequence generator 212 computes a matched 
noise sequence (i.e., a sequence with the same parameters as an observed signal block 
using the parameters estimated by the parameter estimation module 206. In effect, the noise 
sequence generator 212 generates synthetic noise e using parameters of the original sample 
of noise input to the parameter estimation module 206. The synthetic noise sequence is a 
sequence of alpha-stable random variables of the same length as the original sequence 
x ( noise input to the parameter estimation module 206 (i.e., a matched noise sequence made 
up of random numbers having an alpha-stable distribution). In one embodiment, the matched 
noise sequence is generated using an alpha-stable random number generator, which is 
known in the art as disclosed by J. M. Chambers, C. L. Mallows, and B. W. Stuck, in "A 
Method For Simulating Stable Random Variables," Journal of the American Statistical 
Association, Vol. 71, No. 354, pp. 340-344, June 1976, and hereby incorporated herein by 
reference. 

After generating a sequence of alpha-stable variables using the random noise 
generator 212, the differencer 214 subtracts this sequence of synthetic noise e from the 
observed signal block x t _j, thereby converting skewed noise into symmetric noise. The 
resulting signal block xj-I output from the differencer 214 is a modified signal block composed 
of a data component and centro-symmetrized (i.e., deskewed and centralized) noise 
component. In effect, subtracting e from x t -l results in the addition of random noise to the 
observed signal blocker/./, thereby making noise the resulting signal block symmetric. 
The modified signal block .v_/./ is then used by the non-linear prediction filter to estimate a 
cleaned signal block v t . Advantageously, the random noise sequence generator 212 and the 
differencer 214 provide an apparatus for centro-symmetrizing impulse noise in an observed 
signal x t -l so that the l p -norm minimization technique for estimating the parameters of the 
Volterra filter is unbiased (at least when a zero th order term is included and when no self-terms 
are included in this Volterra filter i.e., terms of the form 6 iJ( c ijk where any pair of ij.k are 
equal). 
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F. Parameter Estimation Module 

Figure 8 illustrates a flow diagram that sets forth the steps that are performed by the 
parameter estimation module 206 to estimate the parameters of an alpha-stable distribution. 
By way of overview, the steps performed by the parameter estimation module can be 
summarized as follows. Initially, at step 802 a block (or sequence) of observed data signals or 
samples S = {Xfc} = {Xq, ... A/,./} is received at the parameter estimation module 206. The 
samples are obtained by observing a signal block when the switches 118 are set to operating 
position B. 

At step 804, a determination is made whether to transform the observed data received 
at step 802. Depending on the outcome of the determination made at step 804, one or more 
transformations are performed on the observed data to obtain deskewed (i.e., symmetric) or 
centralized alpha-stable random variables at step 806. Once the transformation of the 
observed data is complete, moments of the alpha-stable distribution are computed at step 
808. Using the computed moments, estimates of the parameters a, p, y, and 5 of an alpha 
stable distribution are computed at step 810. 

Step 804 is repeated depending upon whether all parameters were estimated at step 
812. Once all parameters of the alpha-stable distribution have been computed, the 
parameters are output to the signal estimation module at step 814. it will be appreciated by 
those skilled in the art that the method set forth in Figure 8 need not be used to compute 
every parameter, but instead can be used to estimate a subset of the four parameters a, p, y, 
and 5 of an alpha-stable distribution. 

F.1 Transformations 

More specifically at step 804, a decision is made whether to perform one or more 
transformations on the sequences of data signals Xk- Those transformations selected to be 
performed at step 804 are computed at step 806. The purpose of performing a transformation 
is to eliminate one or more of the parameters in the distribution, thereby minimizing the 
number of variables that are being solved for at any one time. The transformations presented 
below in Tables 1-4 are used to generate, for example, sequences with 5=0 or p=0, or with 
sequences with both 5=0 and p=0 (except when ct=1). Advantageously, by using such 
sequences, methods that can be applied to symmetric variates can be applied to skewed 
variates. In addition, skew-estimation methods for centered variates can be applied to non- 
centered variates, with a loss of some sample size. 

The transformations that can be performed at step 804 include a centro- 
symmetrization transformation Xk cs , a symmetrization transformation xfi l and a 
centralization transformation Xj£. Another available transformation at step 804 is a relocated 
or approximately centralized transformation Xjfi t which requires an estimate of the location 
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parameter 6. Each of these transformations are set forth in Tables 1-4, respectively. More 
specifically, each of the Tables 1-4 set forth a particular transformation that takes the 
weighted sums of the sequences of noise signals Xk (i.e., sequence of stable variates). 

The resulting transformed sequence of noise signals are defined in Tables 1-4 in the 
form of the parameters of the alpha stable distribution S a (dispersion parameter y, symmetry 
parameter p, location parameter 5) for some value of the characteristic exponent a (e.g., a = 
1.5). Which one or ones of these four transformations are performed at step 806 depends on 
the particular variable of the alpha-stable distribution being solved at step 810. 



Table 1: Centro-Symmetrization Transformation 



Transformation 


Resulting Parameters 


x£ s = X lk - X 2 k-\ 


S a (2y, 0, 0) 



Table 2: Symmetrization Transformation 



Transformation 


Resulting Parameters 


xf =X U + X lk _ x -2 Va X u _2 


S a (4y, 0, [2-2 Ua ]S) 



Table 3: Centralization Transformation 



Transformation 


Resulting Parameters 






So([2+2a]X) ^ /?,0) 




Table 4: Relocation Transformation 






Transformation 


Resulting Parameters 






x k = X k - S 







F.2 Computing Moments of Alpha-Stable Distributions 

After transforming the observed noise signals at step 806 if necessary, moments for 
the alpha-stable distribution are estimated at step 809. Estimating a moment of an alpha- 
stable distribution involves evaluating the equations set forth in Tables 5-10 with n samples 
(where L= n samples in Figure 2) of the transformed signals. It will be appreciated by those 
skilled in the art that other choices of sample length may be employed to provide a better 
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trade-off between computation or sampling time and variance of the parameter estimates. 
Specifically, set forth in Tables 5-10 are formulas for computing up to six different classes of 
moments. The different classes of moments include absolute fractional lower order moments 
(FLOM), signed FLOM, signed logarithmic moments, absolute logarithmic moments, extreme 
value moments, and empirical characteristic function moments each of which is set forth in 
Tables 5-10 respectively. In these tables, random variables are denoted by capital letters 
(e.g., example X ). 

In Tables 5 and 6, the p th order moment should be chosen on the basis of a lower 
bound on the possible value of the parameter alpha. If this lower bound is cc min , then a value of 
P=<Xmir/4 is a good choice of p. The value of p should not be chosen to be too large, since if p 
is greater than a/2 then the variance of the FLOM estimator is infinite, and the variance of the 
alpha estimate is therefore large. If p is too small then the absolute FLOM will be dose to one, 
and the variance of the alpha estimate again becomes large. 



Table 5: Absolute FLOM 



Moment 


Moment Estimate 




1 " 

A p =~Y} x k\ P . where p 
n k=\ 

is the order of the 
moment 



Table 6: Signed FLOM 



Moment 


Moment Estimate 


S p = EX <P> 


S p =±£sign(X k )\X k \ p , 
n k=\ 

where p is the order of the 
moment. 




Table 7: Signed Logarithmic 



Moment 


Moment Estimate 


A=.£sign(X)log|X| 


A=-Ssign(^)log|^| 



Table 8: Absolute Logarithmic 



Moment 


Moment Estimate 


L x =£log|X| 

L 2 =£(log|X|-i,) 2 

L 3 =£(log|X|-I,) 3 


4=- L T (XH^|-A) 2 

n 1 k=\ 
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Table 9: Extreme Value 



Moment 


Moment Estimate 


Kj = E max log X 
K, = £max log(-X) 

Y 2 = £((maxlogX)- Y,) 2 
Y 2 = EttmaxlogC-X))^) 2 


Initially, find the maximum (denoted by over 
bar) and minimum (denoted by underscore) 
logarithm in each length r block of the data 
sample, thus: 

K k = max{log Jf r(k _ 1)+l ,log AV (k _ 1)+2 ,--* 

lo E ^r(k-1) + r-l) 

K k = max{log -^ r(k _ 1)+l , log -A; (k _ 1)+2 ,-.- 

lo g -^r(k-1) + r-I> 

Subsequently, compute the mean and 
variance of these Kto obtain the estimates 
of the Y's, thus: 

. n f r 

^ = — Y K k 

. n/r 

Y, = -K) 2 



Table 10: Empirical Characteristic Function 



Moment 


Moment Estimate 


<p(p) = £e^ x 


<p(p) = —^ e JpXk .where 

p is the order of the 
moment. 



F.3 Estimating Parameters Using The Computed Moments 

Using the moments computed using Tables 5-10, the parameters of an alpha-stable 
distribution a, p, y, and 5 are computed using the formulas given in Tables 11-14. Each of the 



Tables 11-14 have an "ID" column, a "condition" column, and an "estimator" column. The "ID" 
column identifies different estimators for the same parameter. The "condition" column defines 
when a particular moment computed at step 808 may be applied. For some of the estimates of 
the parameters, there is included a lower bound on the alpha estimate (i.e., a mjn ). For these 
cases, the a min prevents numerical problems from arising and improves the performance of 
the estimators in situations where such a bound is available. It has been found that a good 
estimate of ot min for signal transmission systems is one (i.e., cc min =1). It will be appreciated that 
dependent on which of the transformations from Tables 1-3 are applied prior to application of 
these estimators, it wiil be necessary to re-transform the estimates obtained for the 
transformed sample back to the parameter values for the original sample. 

Some of the estimators in the Tables 11-14 include a superscript X or Y on the 
moment as in estimators a 2 and <J 3 set forth in Tables 11B and 11C. The presence of a 
superscript X or Y means that the noise samples (e.g., signal block x t ) are partitioned into two 
parts U and V, with each part containing data samples U 1( U 2 , U 3 , ... and V 1( V 2t V 3 , 
respectively. The moment with superscripts is computed for the summed samples as: 

X^IVV,, X 2 =U 2 +V 2 , X 3 =U 3 +V 3 

while that for superscript Y is computed for the concatenated samples as: 

Y^U,, Y 2 =V 1f Y 3 =U 2 , Y 4 =V 2 , Y 5 =U 3 , Y 6 =V 3 

In addition, some of the estimators of alpha in the Tables 11A-11E include an 
auxiliary variable Z. The auxiliary variable Z is used to denote some intermediate function of 
certain moments to simplify the exposition. Also in the estimator d x set forth in Table 1 1 A, the 
arcsinc function is used. By definition, the arcsinc function is the inverse of the sine function 
(i.e., if y = sinc(x) = sin(x)/x and if 0 < x < ti, then x = arcsinc(y)). In the estimator S x in Table 
1 3, a sample's f-fractile is computed. A sample's f-fractile is the point x for which a fraction f of 
the sample lies below x. For example, the lower quartile of the data is the 0.25-fractile and the 
median is the 0.5-fractile. 

In the estimator S 2 in Table 13, the p% truncated sample mean is computed. The p% 
truncated sample mean is the mean of all samples except those larger than the (p/2)% largest 
and smaller than the (p/2)% smallest samples. For example, given a sorted list of one- 
hundred samples, the p% truncated sample mean is computed by truncating p/2 of the largest 
and p/2 of the smallest samples in the sorted list of samples. 
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Table 11 A: Estimate of Alpha 



ID 


Condition 


Estimator 




S = 0, 
/? = <>, 

< «min 
2 


. tan(f) 
2 = 2 

If Z > sinc(7rp/a min ), then a = a min , 
else if Z < sinc(7rp/2), then a = 2 
otherwise a = ;zp(arcsinc(Z))~\ 



Table 11B: Estimate of Alpha 



ID 


Condition 


Estimator 




<? = 0, 

2 


Z = Iog^ -log^ 

If Z<-^ , then a =a min , 

else if Z > ^ ° , then a = 2 

2 

. . pIog2 

otherwise a = ==— . 

Z 


Table 11C: Estimate of Alpha 


ID 


Condition 


Estimator 


«3 


£ = 0 


Z = Zf -if 

- los^ 
If Z then a =a min , 

^min 

else if Z > — — , then a = 2 
2 

loo 9 

otherwise a = — ~. 

Z 
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Table 11D: Estimate of Alpha 



ID 


Condition 


Estimator 


<*4 


£ = 0 


- r L r" 3 

Z = 1 — - , where 

r i 

^ = 1.2020569- = ArHx) . 

l f ^<^min> then ^=^min. 

else if Z > 2, then a - 2 
otherwise a = Z. 



Table 11E: Estimate of Alpha 



ID 


Condition 


Estimator 


^5 


<5 = 0 


Z = — =r + — 
2V6U2 I 2 J 

else if Z > 2, then a = 2 
otherwise a = Z. 



Table 12: Estimates of Beta 



ID 


Condition 


Estimate 




J = 0, estimate 
of a available 


~ _ tan(S 0 a;z/2) 
tan(cuz/2) 


h 


(5 = 0, estimate 
of a available 


^_ tantfA/ZOcwz??)) 
tan(a;z/2) 


k 


5 = 0, estimate 
of a available 


1 . 

expta^i-I,)) 
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Table 13: Estimates of Gamma 



ID 


Condition 


Estimate 


Y\ 


5 = 0, estimate 
of a,J3 available, 
p < a/2 


r-{ r(1 "^ )cos( ^ /2) flcos^i, 

\T(\- pla)cos(p9Ia) ) 
where 6 - arctan(/? tan(a^ / 2)). 


h 


5 = 0, estimate 
of a, /? available 


y - exp(al^ + ^(1 -a))\ cos 6 1, 
where 0 = arctan(/? tan(cr;r / 2)) 
and where y/ is as defined in Table 1 ID . 



Table 14: Estimates of Delta 



ID 


Condition 


Estimate 




Estimate of a,/? 
available 


S = the / -fractile of the sample, where 

1 0 

f = , 0 = arctan(/?tan(a7r/2)). 

2 na 


$2 


J3 = 0 


8 = the p% truncated sample mean 
(this also includes the median). 



F.4 Origins Of The Parameter Estimators For Stable Distributions 

Sections F.4.1-F4.4 describe the principles used to derive the equations set forth in 
Tables 1-14. 

F.4.1 FLOM Estimators 

The estimators based on the fractional lower order moments (FLOM) are all 
rearrangements of the formula in Theorem 1. 



Theorem 1: If X is a stable random variable with parameters a,j3,y, and with £ = 0 



then: 



£[(sign(A')) A W] = 



cos# 



cos(^-if) 



r(i - P ) 

with p e (-2,-1) u (-l,a), p * 1, for k = 1, 
or with p e (-!,«), p * 1 for k = 0, 



cos(^) ' 
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where tan# = /?tan^. 

The proof of Theorem 1 is disclosed by V. M. Zolotarev, in "One-dimensional Stable 
Distributions", Providence, Rl: AMS, 1984. 

F.4.2 Logarithmic Estimators 

The estimators based on logarithmic moments are the consequence of differentiating 
the formula of Theorem 1 and rearranging the formulae obtained by applying the following 
result: 

Lemma 2: Assuming the necessary derivatives exist for a random variable X, 

£[(log|^|)*] = lim £[|AT]for* = l,2,-... 
/>-»<) dp k 

The proof of Lemma 2 is arrived at by differentiating the moment generating function for the 
logarithmic process. 

F.4.3 Extreme Value Estimators 

Extreme value estimators are parameter estimators for the Frechet distribution which 
the tails of the stable pdf obey, which is given by following theorem: 

Theorem 3: The tails of a stable pdf are asymptotically described by: 

lim Pr(X > A)- C(a)(l + P)X~ a 
lim Pr(X < -A.) - C(a)(\ - J3)(-Ay a 
for a suitable function C(a). 

The proof of Theorem 3 is disclosed by G. Samorodnitsky and M. S. Taqqu, in Stable Non- 
Gaussian Random Processes: Stochastic Models with Infinite Variance, Chapman & Hall, 
New York, 1994. 

F.4.4 Weighted Empirical Characteristic Function Estimators 

The empirical characteristic function estimator has been described by S. Kogon and 
D. B. Williams, in "On The Characterization Of Impulsive Noise With Alpha-Stable 
Distributions Using Fourier Techniques," Asilomar Conf. on Signals, Systems, and 
Computers, 1995. The weighted version of this estimator may be derived by: 

1) Taylor expanding the residuals in the classical characteristic function method up 
to 1 st order in the moment estimates. 

2) Approximating the covariance matrix of the residuals using this Taylor expansion. 



-24- 



3) Computing the maximum likelihood estimate of the parameters assuming the 
residuals have a normal distribution with covariance described by the approximated 
covariance matrix. 

G. Estimation Examples 

Figure 9A is a flow diagram which sets forth one combination of steps for estimating 
the parameters a, p, y, and 5 of an alpha-stable distribution. More specifically, Figure 9A sets 
forth an example of a particular sequence in which the steps in Figure 8 can be performed. It 
will be appreciated by those skilled in the art that the flow diagram sets forth only one of many 
different possible sequences in which the estimators in Tables 11-14 can be applied, as 
evidenced by the condition column in the Tables. 

Initially at step 902, a data sample S (e.g., signal block x t .j) is observed with the 
switch 118 in operating position B. At step 904, the centro-symmetrization transform set forth 
in Table 3 is applied to the data sample S to obtain a transformed data sample C. A 
determination is made at step 906 whether a lower bound (i.e., a min ) on alpha is known. !n one 
embodiment, the lower bound on alpha is assumed to equal one - this is an appropriate 
choice for most communication systems, if there exists such a lower bound on alpha, then an 
estimate for the alpha parameter is computed at step 908 by applying the alpha estimator a 2 
to the data sample C; otherwise, the alpha parameter is computed at step 910 by applying the 
alpha estimator a 3 to the data sample C. The alpha estimators a 2 and a 3 are defined above 
in Tables 1 1 A and 1 1 B, respectively. 

To estimate the parameter 5 for the data sample S, steps 912 and 914 are performed. 
At step 912, the symmetrization transformation set forth in Table 2 is applied to the data 
sample S to obtain a transformed data sample T, At step 914, the parameter 5 is estimated as 

the median of the transformed data sample T using the delta estimator S 2 , set forth in Table 
14. This estimate is divided by (2-2 1/Gr )to obtain a delta estimate appropriate for the 
untransformed sample. Subsequently at step 916, the data sample S is relocated using the 
estimate of delta computed at step 914 to obtain a transformed data sample R. At step 918, 
the parameter beta is estimated by applying the beta estimator p x set forth in Table 12 to the 
transformed data sample R. In addition at step 920, the parameter gamma is estimated by 
applying the gamma estimator y x set forth in Table 13 to the transformed data sample R. At 
step 922, the estimated parameters for the alpha-stable distribution are output to the signal 
estimation module at step 922. 

In another embodiment, Figure 9B sets forth another combination of steps 950-961 
that can be performed to estimate the parameters a, p, y, and 5 of an alpha-stable distribution. 
More specifically, the estimator of the parameters is a weighted empirical characteristic 
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function estimator (see Table 10). It is known to perform an empirical characteristic function 
method without performing the steps 955, 956, 960, and 961. Advantageously, the additional 
steps 955, 956, 960, and 961 greatly reduce the variance of the estimates of the parameters. 
It will be appreciated by those skilled in the art that it may be advantageous to iterate through 
steps 955, 956, 960, and 961 more than once to yield better estimates. 

The embodiment shown in Figure 9B first observes a data sample S (e.g., signal 
block .v,_/) with the switch 118 in operating position B. Since the characteristic function is the 
Fourier transform of the pdf of a distribution, it is necessary to select some frequencies (i.e. 
arguments of the characteristic function) to use for estimation. At step 950, the set of 
frequencies [r,,/ 2 ■■•-,/,„] is chosen to be a sequence of positive real numbers. The numbers 
are chosen to be positive in order to simplify the presentation of subsequent steps of the 
estimation procedure. However, it is important that the numbers be unique and non-zero. A 
good choice for these numbers has been found to be [0.05, 0.10, 0.15, 0.90, 0.95, 1.0]. 

At step 952, the centro-symmetrization transform set forth in Table 3 is applied to the 
data sample S to obtain a transformed data sample Y. At step 953, the empirical characteristic 
function is estimated at each of the frequencies selected at step 950, using the formula given 
in Table 10. The logarithm of the logarithm of the characteristic function estimate at frequency 
t k is computed and assigned to a variable y/ k . From the formula for the characteristic 
function of an alpha-stable random variable (i.e., equation (1)), it can be seen that such a 
double logarithm has a linear dependence on the characteristic exponent a and on the 
logarithm of the dispersion log/. Hence a linear regression is used to estimate these 
parameters. 

Since the residuals in this regression are correlated, good estimates are not expected 
unless a weighting matrix is employed to decorrelate them. However,' the extent of the 
correlation depends on the values of the characteristic exponent and dispersion parameters, 
which are what is being estimated. Therefore, an iterative solution procedure is employed in 
which the weighting matrix and the parameters are alternately estimated. The solution 
procedure is initialized at step 951 by assuming that the weighting matrix is the identity matrix. 
New parameter estimates are obtained at step 954. Using these parameters a new weighting 
matrix is determined at step 955. At step 956, a more accurate set of parameter estimates is 
produced. It is possible to iterate this procedure a number of times. However, it has been 
found that a single iteration (as shown in Figure 9B) usually provides most of the improvement 
in the estimates that can be obtained. 

Next it is necessary to estimate. the skew and location parameters of the distribution. 
This is accomplished at steps 957-961 by making use of the characteristic exponent and 
dispersion estimates obtained at step 956. At step 957, the imaginary parts of the logarithm of 
the empirical characteristic function are computed for the original data (rather than the centro- 
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symmetrized data) using the formula given in Table 10. These quantities are divided by their 
frequency and assigned to the variables co k . From the formula for the characteristic function 
of an alpha-stable random variable (i.e., equation (1)), it can be seen that these quantities 
have a linear dependence on the skew parameter /? and on the location parameter Hence 
step 957 performs a linear regression to estimate these parameters. 

The regression is again performed iteratively, starting from an identity matrix estimate 
of the weight matrix at step 958 and producing an improved estimate of the weight matrix at 
step 960. The formula for the weight matrix has been given in terms of the real and imaginary 
parts of the characteristic function evaluated at the frequencies chosen in step 950 and at the 
sums and differences of these frequencies. 

Finally, after one or more iterations, at step 961 the estimated parameters for the 
alpha-stable distribution are output to the signal estimation module 200. 

H. Alternate Operating Environment 

Figure 10 illustrates an alternate operating environment for performing the present 
invention. The operating environment illustrated in Figure 10 is directed at an image 
processing system (i.e., signal processing system), and more particularly to an image 
processing system for cleaning impulse noise from digitally recorded or artificially generated 
images. Figure 10 shows hardware components 1012 and software components 1010 of the 
digital image processing system operating on a general-purpose computer 1002. 

When operating, the general-purpose computer 1002 receives digital images from an 
imaging device 1004 or an imaging synthesizer 1028. The imaging device or imaging 
synthesizer may be operating local to the general purpose computer 1002 or may be 
operating on a network 1004 such as the Internet, thereby requiring digital images to be 
received through a transmission medium 1005 coupling the network 1004 and network I/O 
1022 of the general purpose computer 1002. 

Also coupled to the general-purpose computer 1002 is a printer 1007 and a display 
1008 for outputting digital images. Additional hardware components 1012 operating in the 
general purpose-computer 1002 include user I/O 1014, memory 1016, CPU 1018, and storage 
1020. In addition, the software components 1010 operating in the general-purpose computer 
include operating system software 1024, filter switch 118, pure noise extractor 1026, and 
noise suppression module 114. 

Figure 11 illustrates a process for cleaning impulse noise from digital images 
recorded by imaging device 1006 or formulated by imaging synthesizer 1028 in accordance 
with the present invention. In the event the imaging device 1006 is operating, an image of an 
original scene 1 102 is recorded with an imaging device such as a scanner, a digital camera, a 
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camera with a frame grabber, or the like. Typically, an image captured by an imaging device 
includes noise inherent in the image signal sampled using the imaging device. 

One source of impulse noise corrupting noisy digital image 1104 is the transmission 
medium 1005. Noise that degrades the quality of sampled image signals can either be signal 
dependent noise or additive noise. It is assumed for the purposes of this invention that 
impulse noise that corrupts image data is additive, similar to impulse noise corrupting data 
signals transmitted over analog data channel 110 (shown in Figure 1). 

The filter switch 118, as set forth above, has two operating positions. The operating 
position A is the normal operating mode of the noise suppression module 114. The elements 
forming the noise suppression module 114 are set forth in Figure 2 and described above. In 
normal operating mode, noisy images are cleaned as described above to produce an estimate 
of a clean image 1 106. In operating position B, the filter switch 118 directs noisy digital image 
1104 to the pure noise extractor 1026. The purpose of the pure noise extractor 1026 is to 
provide the parameter estimation module 206 with an observed signal block that consists 
entirely of impulse noise that is absent of image content. 

The pure noise extractor 1026 is necessary because the impulse noise corrupting an 
image recorded with the imaging device 1006 cannot be measured independent of the data 
signals. That is, although the noise is additive, it cannot be independently measured as shown 
in Figure 1 using null modem 116. As set forth above, estimating the parameters of an alpha 
stable distribution of the imaging device 1006 requires an observed signal block consisting of 
impulse noise to be input to parameter estimation module 206. 

Figure 12 is a flow diagram that sets forth the steps performed by the pure noise 
extractor 1026 to produce an observed signal block that consists of impulse noise and no 
image content. Initially at step 1202, three images I 1( I 2 , and I 3 of the same original scene 
1102 are recorded using imaging device 1006 and transmitted, if necessary, through 
transmission medium 1005 to pure noise extractor 1026. At step 1204, the difference between 
two of the images recorded at step 1026 is computed (e.g., I, - I 2 ) to define a centro- 
symmetrized difference image. This operation has the effect of canceling the data component 
and performing the centro-symmetrization transformation set forth in Table 1 on the noise 
component. 

At step 1206, an estimate of the characteristic exponent a is obtained by applying one 
of the alpha estimators set forth in Tables 11A-11E to the centro-symmethzed difference 
image. Subsequently at step 1208, a centralizing transformation, which is set forth in Table 3 t 
is applied .to the three images l n , I 2 , and I 3 to define a centralized difference image I 4 . At step 
1210, the centralized difference image I 4 computed at step 1208 is input to the parameter 
estimation module 206, The parameter estimation module 206 computes the parameters^ an 
.alpha-stable distribution by considering each pixel of the image as an independent sample. 
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Once computed, these parameters are then input to the signal estimation module 200 for 
estimating clean image 1 106. 

In a first alternate embodiment of the pure noise extractor 1026, an image recorded 
with the imaging device 1006 or the like that consists of characters or line segments. A 
segment of the image that has no characters or line segments is isolated. Because the 
isolated area is absent image content, it can be input into the parameter estimation module 
206 to estimate the alpha-stable parameters. 

In a second alternate embodiment of the pure noise extractor, a search is performed 
to identify an area of an image recorded with the imaging device 1006 that is smooth or flat. A 
smooth or flat region is one which has a constant background region or that has small 
changes in gray level or luminance across an area of the recorded image. Properties of a 
such a region in an image can be discovered by moving a window over the recorded image 
and detecting when greater than ninety percent of the gray values lie within plus or minus 
epsilon of some gray-value, where epsilon is a pre-selected threshold value. All the points In 
the discovered region are treated as independent samples of an alpha-stable distribution and 
input to the parameter estimation module 206. 

In an alternate embodiment of the non-linear prediction filter 21 0, the signal estimation 
module 200 shown in detail in Figure 2 is configured to accept signal blocks that are two 
dimensional matrices. In the event a one-dimensional Volterra filter is used to estimate a 
cleaned signal block, images recorded by the imaging device 1006 or formulated by the 
imaging synthesizer 1028 are treated as one-dimensional vectors. To accommodate two- 
dimensions in this alternate embodiment, the non-linear prediction filter 210 is modified as 
described below. 

An example of a two-dimensional (i.e., 2-D) non-linear prediction filter 210 is a 2-D 
Volterra system that can be described as: 

v( m 9 ») = X Z a(A j)x{m - n - ■/) + X S S X! * (/ ' *' l)x(m ~ n - jw™ - *> * - o ( 3 ) 

' ./ i j k f 

+ Z ••• ^c(i,j\kj y s,r)x(n7 - i,n- j)x(m - k,n - l)x{m - r,n - s) . 

/ r 

More details of this Volterra model are described in "A Computational Method For The Design 
Of 2-D Nonlinear Volterra Models," by G.F.Ramponi, G.LSicuranza, W. Ukovich, IEEE Trans. 
On Circuits and Systems, Vol. 35, No. 9, September 1988, pp. 1095-1102, which is hereby 
incorporated by reference. 

The 2-D Volterra model set forth above extends up to third order non-linearity. 
Extending the 2-D Volterra model to a fourth order non-linearity provides little improvement in 
noise suppression but is much more computationally intensive. The summations in the 2-D 
Volterra model apply to a neighborhood of the pixel under consideration. For simplicity, only 
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the nine pixels that make up a 3x3 square centered at the pixel being considered are included 
in the sum. However, alternative neighborhood structures can also be applied. In addition, it 
will be appreciated by those skilled in the art that techniques are available for eliminating 
insignificant coefficients in the summations as described by K. C. Nisbet, B. Mulgrew, and S. 
McLaughlin, in "Reduced State Methods In Nonlinear Prediction," Signal Processing, Vol. 48, 
pp. 37-49, 1996. R. D. Nowak and B. D. van Veen, "Reduced Parameter Volterra Filters," 
Proceedings of ICASSP-95^ Vol. 3, pp. 1569-1572, 1995. 

Also, in this alternate embodiment, the matrices in the coefficient optimization module 
208 are constructed by placing every coefficient in the summation of equation (3) in a vector. 
The data terms *(.). *(>(.) and *(.)x(.)x(.) are placed in the vector according to the scheme of 
the 1-D embodiment. This produces a matrix equation for the 2-D embodiment identical in 
form to the 1-D embodiment. The only difference between the 1-D and 2-D embodiments is 
that the entries of the vector of coefficients are defined by the above-described neighboring 
structure. Once complete, the coefficient optimization module 208 is run as described above 
for the 1-D embodiment. 

I. Summary 

It will be appreciated that the present invention may be readily implemented in 
software using software development environments that provide portable source code that 
can be used on a variety of hardware platforms. Alternatively, the disclosed system may be 
implemented partially or fully in hardware using standard logic circuits. Whether software or 
hardware is used to implement the system varies depending on the speed and efficiency 
requirements of the system and also the particular function and the particular software or 
hardware systems and the particular microprocessor or microcomputer systems being 
utilized. 

The invention has been described with reference to a particular embodiment. 
Modifications and alterations will occur to others upon reading and understanding this 
specification taken together with the drawings. The embodiments are but examples, and 
various alternatives, modifications, variations or improvements may be made by those skilled 
in the art from this teaching which are intended to be encompassed by the following claims. 
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CLAIMS: 

1. A method for reducing impulse noise in a signal processing system, comprising the 

steps of: 

estimating parameters of an alpha-stable distribution to model impulse noise that 
corrupts data signals input into a transmission medium of the signal processing system; 

sampling signals from the transmission medium; said sampling step storing the 
sampled signals in a memory; the sampled signals having a noise component and a data 
component; and 

computing with a prediction filter an estimate of the data components of the sampled 
signals using the estimated parameters of the alpha-stable distribution. 

2. The method according to claim 1, wherein said sampling step is performed from a 
digital subscriber line (DSL). 

3. The method according to claim 2, further comprising the steps of: 

recording samples of impulse noise signals ,v 2 , ... , x n } transmitted over the 
transmission medium of the signal processing system; the samples of impulse noise signals 
recorded by said recording step having no data component; and 

using the recorded samples of impulse noise to perform said estimating step. 

4. The method according to claim 1, wherein said sampling step is performed using 
an imaging device. 

5. The method according to claim 4, further comprising the steps of: 

recording a first, a second, and a third image with the imaging device; 

computing a difference between the first and the second image to define a centro- 
symmetrized difference image; 

estimating a characteristic exponent of an alpha-stable distribution using the centro- 
symmetrized difference image; and 

computing a sample of impulse noise signals by applying a centralizing transformation 
to the first, the second, and the third images. 

6. The method according to claim 1, further comprising the steps of: 

generating a synthetic noise signal with a random number generator that uses the 
parameters of the alpha-stable distribution as input parameters to the random number 
generator; and 
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subtracting the synthetic noise signal from the sampled signal to produce a modified 
sampled signal; the modified sampled signal having a noise component with a symmetric 
distribution. 

7. The method according to claim 1, further comprising the steps of: 

estimating a characteristic exponent of the alpha-stable distribution that models 
impulse noise corrupting the sampled signal; and 

optimizing model coefficients of the prediction filter using the estimated characteristic 
exponent of the alpha-stable distribution. 

8. The method according to claim 7, wherein said optimizing step minimizes a p th - 
power error criterion to optimize the model coefficients of the prediction filter. 

9. The method according to claim 8, wherein said optimizing step defines the p th - 
power error criterion using the characteristic exponent of the alpha-stable distribution. 
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10. The method according to claim 8, wherein said optimizing step minimizes the p lh - 
power error criterion by performing the following computation iteratively: 

where C are model coefficients, W is a diagonal weight matrix, and x t is an 
observed signal block of L samples, where: 
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x 2 [txL + N] x[txL + N]x[txL + N-l] ••• 
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x'[rxL + \] 



-(3) 



0 
0 

r[/xZ] 
x 2 [txL + L-N] 

0 
0 



.x'[/xI + 1].y[/xZ] 
A- 3 [rxZ, + 7V] x 2 [txL + N]x[txL + N~\] ,y 3 [/xL] 
.Y 3 [r x L + L-\] x 2 [t xL + L- l]x[t x L + L - 2] • ,v 3 [r x Z, + L - A'] 

11. The method according to claim 8, wherein said step of estimating the 
characteristic exponent comprises the steps of: 

computing a moment of the alpha-stable distribution; and 

estimating the characteristic exponent using the computed moment. 
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12. The method according to claim 11, wherein said step of computing the moment of 
the alpha-stable distribution comprises the step of computing an absolute fractional lower 
order moment A p according to the following equation: 

A p =- n i \X k \" , where p is the order of the moment and x k is one sample from a 
sequence of noise samples {X u X 2 X n ). 

13. The method according to claim 12. wherein said step of estimating the 
characteristic exponent comprises the step of estimating the characteristic exponent by 
computing z = log a x p -log A r p using the computed absolute fractional lower order moment j, 
by performing the steps of: 

partitioning the sequence of noise samples into two parts U and V, with each part 
containing data samples U, U 2 , U 3 , ... and V,, V 2 , V 3 respectively; 

computing the moment ^by summing noise samples in the sequence of noise 
samples as: 

X^U^V,, X 2 =U 2 +V 2 , X 3 =U 3 +V 3 , and 

computing the moment i'by concatenating noise samples in the sequence of noise 

samples as: 

Yi=U,. Y 2 =V„ Y 3 =U 2 , Y 4 =V 2 , Y 5 =U 3 , Y 6 =V 3 

14. The method according to claim 13, wherein said estimating step estimates the 
characteristic exponent to be equal to -£i2ii when i < £jog2 and f > plog2 _ whefe jg ^ 

order of the moment. 

15. The method according to claim 1 1 , wherein said step of computing the moment of 
the alpha-stable distribution comprises the step of computing an absolute logarithmic moment 

using the following equation: 

L x = -g iog|A-,| , where X k is one sample from a sequence of noise samples {X u X 2 , 

U 

16. The method according to claim 15, wherein' said estimating step estimates the 
characteristic exponent by computing Z = L? -L\ using the computed absolute logarithmic 
moment Z, , by performing the steps of: 
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partitioning the sequence of noise samples into two parts U and V, with each part 
containing data samples U 1( U 2( U 3( ... and V 1( V 2( V 3 , ... ? respectively; 

computing the moment Lf by summing noise samples in the sequence of noise 
samples as: 

X^U^V,, X 2 =U 2 +V 2 , X 3 =U 3 +V 3l and 

computing the moment Zf by concatenating noise samples in the sequence of noise 
samples as: 

Y 1= U 1( Y 2 =V 1( Y 3 =U 2l Y 4 =V 2 , Y 5 =U 3 , Y 6 =V 3 

17. The method according to claim 16, wherein said estimating step estimates the 
characteristic exponent to be equal to i^f! when Z<^^ and Z>^^- 

18. The method according to claim 16, wherein said computing step computes a 
signed logarithmic moment using the following equation: 

A = -Zsign(^)log|ATJ. 
n k=\ ' ' 

19. The method according to claim 18, further comprising the step of estimating a 

symmetry parameter of an alpha-stable distribution by computing tan « A/x i) ^ 

tan(ar;z/2) 

20. The method according to claim 1, wherein said computing step extrapolates with 
the prediction filter an estimate of the data components of the sampled signals using the 
estimated parameters of the alpha-stable distribution. 

21. The method according to claim 1, wherein said computing step interpolates with 
the prediction filter an estimate of the data components of the sampled signals using the 
estimated parameters of the alpha-stable distribution. 

22. The method according to claim 1, further comprising the step of computing a 
moment of the alpha-stable distribution. 

23. The method according to claim 22, wherein said step of computing a moment 
comprises the step of computing one of: 

an absolute fractional lower order moment A = — £ \x k \ p , 

rtk=\ 

a signed fractional lower order moment S = - £ sign(A^) \X k \ p , 
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a signed logarithmic moment A = 1 £ sign^^logprJ 



1 

a first absolute logarithmic moment Z, =— ^ log|A^ 



a second absolute logarithmic moment L 2 =— I— (^T i 0 g|^|- Zi) 2 



a third absolute logarithmic moment Z 3 =— — logpr^-i,) 3 , 

an empirical characteristic function moment <p(p) = 1 V e >^ 

a first upper extreme value moment Y x =— — Vtf* 

a first lower extreme value moment Y x = — L_ VV 

a second upper extreme value moment F, = ! Vf7j _vV 

(/i/r)-l£ l i; ' 

1 

a second lower extreme value moment Y 7 = Vr/r -9 ^ 2 and 

( w / r )_!f,^ il' ' dna 

/r=l 

where A> is one sample from a sequence of noise samples {X 1t X 2} ... ,x n }, 
p is the order of the moment, 

K k =max{logJr r(k . 1)+I ,iog^ r(k . 1)+2 ,..- log J r( , 1)+r , } , 
.K k =max{log -^ r(k _ 1)+l , log -^ r(k . 1)+2 log -I r(k , )+H } , and 
r is a block length of data. 

24. The method according to claim 23, further comprising the step of computing a 
characteristic exponent a of the alpha-stable distribution. 

25. The method according to claim 24, wherein said step of computing the 
characteristic exponent a of the alpha-stable distribution comprises the step of computing one 

of: 

a first auxiliary variable Z = \o°A* -\ogA Y , 
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a second auxiliary variable Z = Zf -L\ , and 

r L y" 3 

a third auxiliary variable Z = 1 - ^ 7Q2 Q5 6 9 1 and 



,1/2 

rr ( \ 1 

a fourth auxiliary variable Z = - 



26. The method according to claim 24, further comprising the step of computing a 
symmetry parameter p of the alpha-stable distribution. 

27. The method according to claim 26, wherein said step of computing the symmetry 
parameter p of the alpha-stable distribution comprises the step of computing one of: 

- tan((A/Z,)a^/2)) 

a first estimate of the symmetry parameter p = — 1 , and 

tan(a^/2) 

2 

a second estimate of the symmetry parameter p = 1 = . 

exp(a(F, -V { )) 

28. The method according to claim 24, further comprising the step of computing a 
dispersion y of the alpha-stable distribution. 

29. The method according to claim 28, wherein said step of computing the dispersion 
y of the alpha-stable distribution comprises the step of computing 

f r(i- p)cos(ptt/2) Y /p 



r = 



{r(\-p/a) cosipOia) 



| cos# | , where 6 = arctan(/? tan(a^ 12)) 



30. The method according to claim 24, further comprising the step of computing a 
location parameter 5 of the alpha-stable distribution. 

31. The method according to claim 30, wherein said step of computing the location 
parameter 5 of the alpha-stable distribution comprises the step of computing a fractile /of the 

1 0 

sequence of noise, where / = , 6 = arctan(/? tan(<27r / 2)) . 

2 Ka 

32. The method according to claim 1, further comprising the steps of: 
computing a characteristic exponent a of the alpha-stable distribution 
computing a location parameter 5 of the alpha-stable distribution, 
computing a symmetry parameter p of the alpha-stable distribution, and 
computing a dispersion y of the alpha-stable distribution. 
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33. The method according to claim 32, wherein said step of computing the 
characteristic exponent a and said step of computing the dispersion y of the alpha-stable 
distribution comprises the step of computing the following matrix equation: 



/=/».*=«/ j^m,k~m 

X Iog(/ > )log( ^ )C ^ X log( '* )C ;* 

j=fi,.k=m j=m,k=iii 

2>g(/*)C 77t Ya C J" 

j=\.k=\ >l,*=l 

where 



a 



j=m.k=ttt 

Ys^kWjCjk 

y=u=i 

j= m,k= tit 
J=\,k=\ 



{/ l ,/ 2 ,-",/ /#l } is some set of positive real numbers, 

C equals W* 1 where W is an mxm matrix with elements 

and a\y % are existing estimates of the characteristic exponent and dispersion, 
and where y/ k = log(- log(Re £(;*))) for sequence r for* = l,--,/w , and 

Y is a centro-symmetrized version of the sequence of noise samples. 

34. The method according to claim 33, wherein said step of computing the location 
parameter 5 and said step of computing the symmetry parameter p of the alpha-stable 
distribution comprises the step of computing the following matrix equation: 



j=mA=m j- in ,k —ni 

./=U=l j=\,k=\ 
j = m.k- in j= m.k= n I 



where k = y tan^^- | , 



j= ni,k= m 

A: =1 
j=ni.k=m 

>U=1 



C equals W" 1 where W is the mxm matrix with elements 
1 



+ l s „(I s R r + R s l r ) + I s _ r {I s R r -R s I r )-AR s l s R r l r ) 



where 

R s = Re(p(/ f )),/, = [m(<p(t s )lR s+r = Re(q>(t s + f r )) f 
where <p(0 = exp(-/ \t f (1 + ^'tan^) + jtS*) for existing estimates p\8* of the skew 



and location parameters, and 
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co k = -\m\og(<p(t k y)f t k for the sequence of noise samples. 
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35. An apparatus for reducing impulse noise in a signal processing system, 
comprising: 

a parameter estimation module for estimating parameters of an alpha-stable 
distribution; the alpha-stable distribution modeling impulse noise that corrupts data signals 
input into a transmission medium of the signal processing system; 

a memory for accumulating sampled signals output from the transmission medium; 
the sampled signals having a noise component and a data component; and 

a signal estimation module for computing an estimate of the data component of the 
sampled signals output from the transmission medium using the estimated parameters of the 
atpha-stable distribution. 

36. The apparatus according to claim 35, wherein said signal estimation module 
further comprises: 

a prediction filter for estimating, using model coefficients, the data components of the 
sampled signals output from the transmission medium; the estimated data components 
corresponding to an estimation of the data signals input into the transmission medium; and 

a coefficient optimization module for optimizing the model coefficients of the 
prediction filter using one of the estimated parameters of the alpha-stable distribution received 
from said parameter estimation module. 

37. The apparatus according to claim 36, wherein the dependence of the prediction 
filter output on the model coefficients is linear. 

38. The apparatus according to claim 37, wherein said prediction filter is a Volterra 

filter. 

39. The apparatus according to claim 36, wherein the dependence of the prediction 
filter output on the data components of the sampled signals input into said prediction filter is 
non-linear. 

40. The apparatus according to claim 36, wherein said parameter estimation module 
further comprises: 

means for computing moments of the alpha-stable distribution; and 

means for estimating the parameters of the alpha-stable distribution with the 
computed moments. 

41. The apparatus according to claim 40, wherein said parameters estimation module 
further comprises means for transforming the alpha-stable distribution to obtain deskewed 
alpha-stable random variables. * 
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